
namespace Tuna {

template<class T_sfield>
void LES_SSF::calcTurbulentViscosity(T_sfield &nut,
                                     const T_sfield &u,
				     const T_sfield &v,
				     const T_sfield &w,
				     double dx, double dy, double dz)
{

    double dc = dx, dc_dx23 = 1.0, dc_dy23 = 1.0, dc_dz23 = 1.0;
    double uip, uim, ujp, ujm, ukp, ukm;
    double vip, vim, vjp, vjm, vkp, vkm;
    double wip, wim, wjp, wjm, wkp, wkm;
    double dxyz = dx * dy * dz;

    dc = pow(dxyz, 1.0/3.0);
    dc_dx23 = pow(dc / dx, 2.0/3.0);
    dc_dy23 = pow(dc / dy, 2.0/3.0);
    dc_dz23 = pow(dc / dz, 2.0/3.0);

    int bi = nut.lbound(firstDim) + 1, ei = nut.ubound(firstDim) - 1;
    int bj = nut.lbound(secondDim) + 1, ej = nut.ubound(secondDim) - 1;
    int bk = nut.lbound(thirdDim) + 1, ek = nut.ubound(thirdDim) - 1;


    for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j) 
	    for(int k = bk; k <= ek; ++k)
	    {
		uip = (u(i,j,k) - u(i+1,j,k));
		uim = (u(i,j,k) - u(i-1,j,k));
		ujp = (u(i,j,k) - u(i,j+1,k));
		ujm = (u(i,j,k) - u(i,j-1,k));
		ukp = (u(i,j,k) - u(i,j,k+1));
		ukm = (u(i,j,k) - u(i,j,k+1));
		
		vip = (v(i,j,k) - v(i+1,j,k));
		vim = (v(i,j,k) - v(i-1,j,k));
		vjp = (v(i,j,k) - v(i,j+1,k));
		vjm = (v(i,j,k) - v(i,j-1,k));
		vkp = (v(i,j,k) - v(i,j,k+1));
		vkm = (v(i,j,k) - v(i,j,k+1));    
		
		wip = (w(i,j,k) - w(i+1,j,k));
		wim = (w(i,j,k) - w(i-1,j,k));
		wjp = (w(i,j,k) - w(i,j+1,k));
		wjm = (w(i,j,k) - w(i,j-1,k));
		wkp = (w(i,j,k) - w(i,j,k+1));
		wkm = (w(i,j,k) - w(i,j,k+1));
		
		nut(i,j,k) = ( uip*uip + vip*vip + wip*wip + 
			       uim*uim + vim*vim + wim*wim ) * dc_dx23 +
		             ( ujp*ujp + vjp*vjp + wjp*wjp + 
			       ujm*ujm + vjm*vjm + wjm*wjm ) * dc_dy23 +
                             ( ukp*ukp + vkp*vkp + wkp*wkp + 
			       ukm*ukm + vkm*vkm + wkm*wkm ) * dc_dz23;

		nut(i,j,k) = 0.104 * dc * sqrt( nut(i,j,k) / 6.0 );
	    }
}

template<class T_sfield>
void LES_SSF::selectTurbulentViscosity(T_sfield &nut,
                                      const T_sfield &v1, 
				      const T_sfield &v2, 
				      const T_sfield &v3,
				      double dx , double dy, double dz)
{
    double cutoff = cos( 20.0 * atan(1.0) / 45.0 );
    double vsix1, vsix2, vsix3, cos_theta, norm_v, norm_vsix, norm, dot;
    
    int bi = nut.lbound(firstDim) + 1, ei = nut.ubound(firstDim) - 1;
    int bj = nut.lbound(secondDim) + 1, ej = nut.ubound(secondDim) - 1;
    int bk = nut.lbound(thirdDim) + 1, ek = nut.ubound(thirdDim) - 1;

    for(int i = bi; i <= ei; ++i)
	for(int j = bj; j <= ej; ++j) 
	    for(int k = bk; k <= ek; ++k)
	    {
		vsix1 = ( v1(i-1,j,k) + v1(i+1,j,k) +
			  v1(i,j-1,k) + v1(i,j+1,k) +
			  v1(i,j,k-1) + v1(i+1,j,k+1) ) / 6.0;
		vsix2 = ( v2(i-1,j,k) + v2(i+1,j,k) +
			  v2(i,j-1,k) + v2(i,j+1,k) +
			  v2(i,j,k-1) + v2(i+1,j,k+1) ) / 6.0;
		vsix3 = ( v3(i-1,j,k) + v3(i+1,j,k) +
			  v3(i,j-1,k) + v3(i,j+1,k) +
			  v3(i,j,k-1) + v3(i+1,j,k+1) ) / 6.0;   

		norm_vsix = vsix1 * vsix1 + vsix2 * vsix2 + vsix3 * vsix3;

		norm_v = v1(i,j,k) * v1(i,j,k) + v2(i,j,k) * v2(i,j,k) +
                         v3(i,j,k) * v3(i,j,k);

		norm = sqrt ( norm_v * norm_vsix );

		dot = v1(i,j,k) * vsix1 + v2(i,j,k) * vsix2 +
		      v3(i,j,k) * vsix3;

		cos_theta = dot / norm ;

// NOTE: here we are comparing the cos(theta) and not the angle theta.
// cos(theta) is a monotonically decresing function, from zero to pi.
		if (cos_theta > cutoff) nut(i,j,k) = 0.0;
	    }
}

} // Tuna namespace
